\section{Derivation of Early Warning Indicators and Composite Index}
\label{app:early-warning}

This appendix provides rigorous derivations of the early warning indicators introduced in Section 3.3, establishes their critical scaling behavior, and develops the optimal composite index construction with statistical validation methodology.

\subsection{Autocorrelation Timescale Divergence}

\textbf{Definition.} The autocorrelation function for order parameter $\Phi_\ell$ is defined as
\begin{equation}
    C_\ell(t) = \frac{\langle \Phi_\ell(t + t_0) \Phi_\ell(t_0) \rangle - \langle \Phi_\ell \rangle^2}{\langle \Phi_\ell^2 \rangle - \langle \Phi_\ell \rangle^2},
    \label{eq:autocorr-def}
\end{equation}
where $\langle \cdot \rangle$ denotes ensemble averaging over initial conditions and stochastic realizations.

\textbf{Exponential decay ansatz.} Far from criticality, relaxation follows exponential decay
\begin{equation}
    C_\ell(t) \sim \exp(-t/\tau_{\text{AC}}),
    \label{eq:exp-decay}
\end{equation}
where the autocorrelation timescale $\tau_{\text{AC}}$ characterizes the slowest decaying mode.

\textbf{Critical scaling from correlation length.} The spatial correlation function near criticality exhibits power-law decay:
\begin{equation}
    G(r) = \langle \phi(x + r) \phi(x) \rangle - \langle \phi \rangle^2 \sim r^{-(d-2+\eta)} \exp(-r/\xi),
\end{equation}
where $\xi$ is the correlation length and $\eta$ is the anomalous dimension. At criticality, the correlation length diverges as
\begin{equation}
    \xi \sim |T - T_c|^{-\nu},
    \label{eq:corr-length-app}
\end{equation}
with critical exponent $\nu$.

\textbf{Dynamic exponent connection.} The dynamic critical exponent $z$ connects spatial and temporal scales through the scaling relation
\begin{equation}
    \tau_{\text{AC}} \sim \xi^z.
    \label{eq:dynamic-scaling-app}
\end{equation}
This follows from the time-dependent Ginzburg-Landau equation:
\begin{equation}
    \frac{\partial \phi}{\partial t} = -\Gamma \frac{\delta \mathcal{F}}{\delta \phi} + \eta(x,t),
\end{equation}
where $\mathcal{F}[\phi]$ is the free energy functional and $\eta$ is thermal noise. Dimensional analysis yields $\tau \sim \Gamma^{-1} \xi^2$ for model A dynamics (non-conserved order parameter), giving $z = 2$.

\textbf{Critical divergence.} Combining \eqref{eq:corr-length-app} and \eqref{eq:dynamic-scaling-app}:
\begin{equation}
    \boxed{\tau_{\text{AC}} \sim |T - T_c|^{-\nu z}}
    \label{eq:tau-divergence-app}
\end{equation}

For the 3D Ising universality class (Class I), $\nu \approx 0.630$ and $z = 2$, yielding $\nu z \approx 1.26$. For mean-field systems (Class III), $\nu = 0.5$ and $z = 2$, giving $\nu z = 1$.

\textbf{Operational measurement.} Estimate $\tau_{\text{AC}}$ from simulation data by fitting exponential decay to $C_\ell(t)$ or computing the integrated autocorrelation time:
\begin{equation}
    \tau_{\text{AC}} = \frac{1}{2} + \sum_{t=1}^{t_{\max}} C_\ell(t).
\end{equation}
Warning threshold: $\tau_{\text{AC}}/\tau_{\text{baseline}} > 5$ where $\tau_{\text{baseline}}$ is measured in the stable regime $T \gg T_c$.

\subsection{Variance Amplification}

\textbf{Fluctuation-dissipation theorem.} In thermal equilibrium, the variance of the order parameter is related to susceptibility $\chi$ by
\begin{equation}
    \text{Var}[\Phi_\ell] = \langle \Phi_\ell^2 \rangle - \langle \Phi_\ell \rangle^2 = k_B T \chi_\ell,
    \label{eq:fdt-app}
\end{equation}
where
\begin{equation}
    \chi_\ell = \frac{\partial \langle \Phi_\ell \rangle}{\partial h_\ell}
\end{equation}
is the response to external field $h_\ell$.

\textbf{Susceptibility divergence.} From linear response theory, the static susceptibility is given by
\begin{equation}
    \chi = \int_V \langle \phi(x) \phi(0) \rangle \, \diff^d x \sim \xi^{2-\eta}.
\end{equation}
Using $\xi \sim |T - T_c|^{-\nu}$, we obtain
\begin{equation}
    \chi \sim |T - T_c|^{-\gamma},
    \label{eq:susc-divergence-app}
\end{equation}
where the susceptibility exponent $\gamma = \nu(2 - \eta)$ follows from scaling relations.

\textbf{Statistical mechanics derivation.} From the Hamiltonian formulation
\begin{equation}
    \mathcal{H}(s) = -\sum_{\langle i,j \rangle} J_{ij} s_i s_j - \sum_i h_i s_i,
\end{equation}
the partition function is $Z = \sum_{\{s\}} e^{-\beta \mathcal{H}}$ and the magnetization
\begin{equation}
    m = \frac{1}{N} \sum_i \langle s_i \rangle = -\frac{1}{\beta N} \frac{\partial \ln Z}{\partial h}.
\end{equation}
The susceptibility becomes
\begin{equation}
    \chi = \frac{\partial m}{\partial h} = \frac{\beta}{N} \left( \langle M^2 \rangle - \langle M \rangle^2 \right),
\end{equation}
where $M = \sum_i s_i$. This establishes $\text{Var}[M] = N k_B T \chi$.

\textbf{Variance scaling.} Combining \eqref{eq:fdt-app} and \eqref{eq:susc-divergence-app}:
\begin{equation}
    \boxed{\text{Var}[\Phi_\ell] \sim |T - T_c|^{-\gamma}}
    \label{eq:var-divergence-app}
\end{equation}

For 3D Ising, $\gamma \approx 1.237$; for mean-field, $\gamma = 1$. Warning threshold: $\text{Var}[\Phi_\ell]/\text{Var}_{\text{baseline}} > 3$.

\subsection{Transfer Entropy Lag}

\textbf{Transfer entropy definition.} The directional information flow from level $\ell$ to level $\ell+1$ is quantified by
\begin{equation}
    TE_{\ell \to \ell+1}(\tau) = \sum_{s_{\ell+1}^{t+1}, s_{\ell+1}^{(k)}, s_\ell^{(k)}} p(s_{\ell+1}^{t+1}, s_{\ell+1}^{(k)}, s_\ell^{(k-\tau)}) \log \frac{p(s_{\ell+1}^{t+1} | s_{\ell+1}^{(k)}, s_\ell^{(k-\tau)})}{p(s_{\ell+1}^{t+1} | s_{\ell+1}^{(k)})},
\end{equation}
where $s^{(k)} = (s^t, s^{t-1}, \ldots, s^{t-k+1})$ denotes the history window and $\tau$ is the time lag.

\textbf{Optimal lag.} Define the characteristic lag $\tau^*$ as the lag maximizing transfer entropy:
\begin{equation}
    \tau^* = \arg\max_\tau TE_{\ell \to \ell+1}(\tau).
\end{equation}
This represents the intrinsic information propagation delay between hierarchical levels.

\textbf{Connection to relaxation modes.} Consider Fourier-space relaxation modes with dispersion relation
\begin{equation}
    \omega(k) = -\Gamma k^2 + i \Delta\omega(k),
\end{equation}
where $\Gamma$ is the kinetic coefficient. The relaxation time at wavevector $k$ is
\begin{equation}
    \tau(k) = \frac{1}{\Gamma k^2} \sim \frac{\xi^2}{\Gamma} \quad \text{for } k \sim \xi^{-1}.
\end{equation}

\textbf{Information propagation speed.} The effective information propagation speed is
\begin{equation}
    v_{\text{info}} \sim \frac{\xi}{\tau(\xi^{-1})} \sim \frac{1}{\Gamma \xi} \sim |T - T_c|^{\nu}.
\end{equation}
Therefore, the characteristic lag for information to traverse spatial scale $\xi$ is
\begin{equation}
    \tau^* \sim \frac{\xi}{v_{\text{info}}} \sim \xi^2 \sim |T - T_c|^{-2\nu}.
\end{equation}

For systems with non-trivial dynamic exponent $z$, the general scaling is
\begin{equation}
    \boxed{\tau^* \sim |T - T_c|^{-\nu z}}
    \label{eq:te-lag-app}
\end{equation}

Warning threshold: $\tau^* > 2\tau_{\text{actuation}}$ where $\tau_{\text{actuation}}$ is the typical actuation delay in the hierarchy.

\subsection{Sobol Sensitivity Explosion}

\textbf{Variance-based sensitivity indices.} For a function $f(X_1, \ldots, X_d)$ with independent inputs, the first-order Sobol index for variable $X_i$ is
\begin{equation}
    S_i = \frac{\text{Var}_{X_i}[\mathbb{E}_{X_{\sim i}}[f | X_i]]}{\text{Var}[f]},
\end{equation}
where $X_{\sim i}$ denotes all variables except $X_i$. This measures the fraction of output variance attributable to variations in $X_i$ alone.

\textbf{Hierarchical context.} For the order parameter $\Phi_\ell$ depending on control parameters $\theta = (\theta_1, \ldots, \theta_d)$:
\begin{equation}
    S_i = \frac{\text{Var}_{\theta_i}[\mathbb{E}_{\theta_{\sim i}}[\Phi_\ell | \theta_i]]}{\text{Var}[\Phi_\ell]}.
\end{equation}

\textbf{Critical amplification.} Near criticality, the order parameter becomes increasingly sensitive to perturbations. The variance in the numerator remains approximately constant (set by parameter ranges), while the denominator diverges as $\text{Var}[\Phi_\ell] \sim |T - T_c|^{-\gamma}$.

Therefore, sensitivity indices scale as
\begin{equation}
    S_i \sim \frac{\text{const}}{\text{Var}[\Phi_\ell]} \sim |T - T_c|^{\gamma}.
    \label{eq:sobol-scaling-prelim-app}
\end{equation}

\textbf{Refined analysis via ANOVA decomposition.} The variance decomposition is
\begin{equation}
    \text{Var}[\Phi_\ell] = \sum_i V_i + \sum_{i<j} V_{ij} + \cdots + V_{1,\ldots,d},
\end{equation}
where $V_i = \text{Var}_{\theta_i}[\mathbb{E}_{\theta_{\sim i}}[\Phi_\ell | \theta_i]]$. Near criticality, coupling between parameters intensifies, causing higher-order interaction terms to grow. However, for parameters directly coupled to the order parameter (e.g., temperature, external fields), first-order terms dominate:
\begin{equation}
    \boxed{S_i \sim |T - T_c|^{\gamma}}
    \label{eq:sobol-divergence-app}
\end{equation}

Warning threshold: $\max_i S_i > 0.7$, indicating that a single parameter accounts for over 70\% of output variance.

\subsection{Composite Index Construction}

\textbf{Individual indicators.} We define normalized versions of the four indicators:
\begin{align}
    I_1 &= \frac{\log(\tau_{\text{AC}}/\tau_0)}{\log(\tau_{\max}/\tau_0)}, \label{eq:I1-app} \\
    I_2 &= \frac{\text{Var}[\Phi_\ell] - \text{Var}_0}{\text{Var}_{\max} - \text{Var}_0}, \label{eq:I2-app} \\
    I_3 &= \frac{\tau^* - \tau_0^*}{\tau_{\max}^* - \tau_0^*}, \label{eq:I3-app} \\
    I_4 &= \frac{S_{\max} - S_0}{S_{\text{crit}} - S_0}, \label{eq:I4-app}
\end{align}
where subscript 0 denotes baseline values far from criticality and subscript max/crit denotes empirically observed maxima or critical thresholds.

\textbf{Composite index formulation.} The early warning system combines these into
\begin{equation}
    \mathcal{W}_{\text{EWS}} = w_1 I_1 + w_2 I_2 + w_3 I_3 + w_4 I_4,
    \label{eq:composite-raw-app}
\end{equation}
equivalent to the form given in equation \eqref{eq:early-warning} after substituting definitions and noting that FDV (fluctuation-dissipation violation) is proportional to $I_2$, and $\Delta TE^{\text{reversal}}$ relates to $I_3$.

\textbf{Principal component analysis for optimal weights.} Given a training dataset with $N$ samples and known proximity to criticality $y_n \in [0,1]$ (distance from $T_c$), we construct the indicator matrix
\begin{equation}
    \mathbf{X} = \begin{pmatrix}
        I_1^{(1)} & I_2^{(1)} & I_3^{(1)} & I_4^{(1)} \\
        \vdots & \vdots & \vdots & \vdots \\
        I_1^{(N)} & I_2^{(N)} & I_3^{(N)} & I_4^{(N)}
    \end{pmatrix}.
\end{equation}

Perform singular value decomposition:
\begin{equation}
    \mathbf{X} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^\top.
\end{equation}

The first principal component $\mathbf{v}_1$ (first column of $\mathbf{V}$) provides the optimal linear combination maximizing variance:
\begin{equation}
    w_i = (\mathbf{v}_1)_i, \quad i = 1,\ldots,4.
    \label{eq:pca-weights-app}
\end{equation}

\textbf{Alternative: ROC-based optimization.} For binary classification (critical vs.\ stable), define true positive rate (TPR) and false positive rate (FPR) for threshold $\theta$:
\begin{align}
    \text{TPR}(\theta) &= \frac{\#\{n : \mathcal{W}_{\text{EWS}}^{(n)} > \theta \text{ and } y_n = 1\}}{\#\{n : y_n = 1\}}, \\
    \text{FPR}(\theta) &= \frac{\#\{n : \mathcal{W}_{\text{EWS}}^{(n)} > \theta \text{ and } y_n = 0\}}{\#\{n : y_n = 0\}}.
\end{align}

The area under the ROC curve (AUC) is
\begin{equation}
    \text{AUC}(\mathbf{w}) = \int_0^1 \text{TPR}(\text{FPR}^{-1}(u)) \, \diff u.
\end{equation}

Optimize weights to maximize AUC:
\begin{equation}
    \mathbf{w}^* = \arg\max_{\mathbf{w}} \text{AUC}(\mathbf{w}) \quad \text{subject to} \quad \|\mathbf{w}\|_1 = 1, \, w_i \geq 0.
    \label{eq:roc-optimization-app}
\end{equation}

This can be solved via gradient-based optimization or grid search.

\textbf{Threshold selection.} Define alert levels based on standard deviations from the mean in the stable regime:
\begin{align}
    \theta_{\text{warning}} &= \mu_{\text{stable}} + 3\sigma_{\text{stable}}, \label{eq:threshold-3sig-app} \\
    \theta_{\text{critical}} &= \mu_{\text{stable}} + 5\sigma_{\text{stable}}, \label{eq:threshold-5sig-app}
\end{align}
where
\begin{align}
    \mu_{\text{stable}} &= \mathbb{E}[\mathcal{W}_{\text{EWS}} \mid T > T_c + \Delta T_{\text{safe}}], \\
    \sigma_{\text{stable}}^2 &= \text{Var}[\mathcal{W}_{\text{EWS}} \mid T > T_c + \Delta T_{\text{safe}}].
\end{align}

The margin $\Delta T_{\text{safe}}$ ensures baseline is measured well away from critical fluctuations (typically $\Delta T_{\text{safe}} \geq 5\sqrt{\text{Var}[T]}$).

\textbf{False positive rate minimization.} Given cost function $C_{\text{FP}}$ for false positives and $C_{\text{FN}}$ for false negatives, choose threshold to minimize expected cost:
\begin{equation}
    \theta^* = \arg\min_\theta \left[ C_{\text{FP}} \cdot \text{FPR}(\theta) + C_{\text{FN}} \cdot (1 - \text{TPR}(\theta)) \right].
    \label{eq:cost-minimization-app}
\end{equation}

For early warning systems, false negatives (missed critical events) typically have much higher cost than false positives (unnecessary interventions), suggesting $C_{\text{FN}}/C_{\text{FP}} \sim 10{-}100$.

\subsection{Statistical Validation Methodology}

\textbf{Bootstrap confidence intervals.} For estimating uncertainty in indicator measurements, perform $B = 10^3$ bootstrap resamples:
\begin{enumerate}
    \item Draw $N$ samples with replacement from the trajectory dataset.
    \item Compute indicators $\{I_1^{(b)}, I_2^{(b)}, I_3^{(b)}, I_4^{(b)}\}$ and composite index $\mathcal{W}_{\text{EWS}}^{(b)}$.
    \item Construct $(1-\alpha)$ confidence interval from percentiles of bootstrap distribution.
\end{enumerate}

The 95\% confidence interval is $[\mathcal{W}_{0.025}, \mathcal{W}_{0.975}]$ where subscripts denote quantiles.

\textbf{Finite-size scaling validation.} To confirm critical scaling \eqref{eq:tau-divergence-app}, \eqref{eq:var-divergence-app}, \eqref{eq:te-lag-app}, \eqref{eq:sobol-divergence-app}, run simulations at system sizes $L \in \{L_1, L_2, \ldots, L_K\}$ and temperatures $T \in \{T_1, \ldots, T_M\}$ near $T_c(L)$.

For each indicator $I$, fit the finite-size scaling ansatz:
\begin{equation}
    I(T, L) = L^{x/\nu} f\left( (T - T_c(L)) L^{1/\nu} \right),
    \label{eq:fss-ansatz-app}
\end{equation}
where $x$ is the scaling dimension of $I$ and $f$ is a universal scaling function. For $\tau_{\text{AC}}$ and $\tau^*$, $x = \nu z$; for Var$[\Phi]$ and $S_i$, $x = \gamma$.

Collapse data by plotting $I L^{-x/\nu}$ versus $(T - T_c) L^{1/\nu}$. Successful collapse confirms critical scaling hypothesis and provides accurate estimates of exponents $\{\nu, z, \gamma\}$.

\textbf{Cross-validation for predictive performance.} Partition dataset into $K = 5$ folds. For fold $k$:
\begin{enumerate}
    \item Train weights $\mathbf{w}^{(k)}$ on folds $\{1,\ldots,K\} \setminus \{k\}$ via \eqref{eq:roc-optimization-app}.
    \item Evaluate AUC and calibration on held-out fold $k$.
    \item Compute Brier score: $\text{BS} = \frac{1}{N_k} \sum_{n \in \text{fold } k} (p_n - y_n)^2$ where $p_n$ is predicted probability of criticality.
\end{enumerate}

Average performance metrics across folds to estimate generalization capability.

\textbf{Threshold calibration procedure.} Given trained composite index, calibrate thresholds using a separate validation dataset:
\begin{enumerate}
    \item For each sample $n$, compute $\mathcal{W}_{\text{EWS}}^{(n)}$ and measure actual distance to criticality $\delta_n = |T_n - T_c|$.
    \item Bin samples by $\mathcal{W}_{\text{EWS}}$ values and compute mean distance to criticality in each bin.
    \item Identify $\theta_{\text{warning}}$ as the value where mean distance drops below $\Delta T_{\text{warning}} \approx 0.1 \times (T_{\max} - T_c)$.
    \item Identify $\theta_{\text{critical}}$ as the value where mean distance drops below $\Delta T_{\text{critical}} \approx 0.05 \times (T_{\max} - T_c)$.
\end{enumerate}

This empirical calibration ensures thresholds correspond to operationally meaningful proximity to criticality.

\textbf{Implementation summary.} The complete early warning system pipeline:
\begin{enumerate}
    \item Collect training data spanning stable and near-critical regimes.
    \item Compute indicators $\{I_1, I_2, I_3, I_4\}$ via \eqref{eq:I1-app}--\eqref{eq:I4-app}.
    \item Optimize weights via PCA \eqref{eq:pca-weights-app} or ROC maximization \eqref{eq:roc-optimization-app}.
    \item Calibrate thresholds using validation data as described above.
    \item During operation, monitor $\mathcal{W}_{\text{EWS}}$ in real-time and trigger interventions when $\mathcal{W}_{\text{EWS}} > \theta_{\text{warning}}$ or $\mathcal{W}_{\text{EWS}} > \theta_{\text{critical}}$.
    \item Periodically retrain weights and recalibrate thresholds as system dynamics evolve.
\end{enumerate}

This methodology ensures statistically rigorous early warning with quantifiable false positive/negative rates, enabling proactive governance before coordination breakdown.
